Ground State of a Homogeneous Bose Gas: A Diffusion Monte Carlo Calculation 
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We use a diffusion Monte Carlo method to calculate the 
lowest energy state of a uniform gas of bosons interacting 
through different model potentials, both strictly repulsive and 
with an attractive well. We explicitly verify that at low den- 
sity the energy per particle follow a universal behavior fixed by 
the gas parameter no? . In the regime of densities typical for 
experiments in trapped Bose-condensed gases the corrections 
to the mean-field energy greatly exceed the differences due to 
the details of the potential. 
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The achievement of Bose- Einstein condensation (BEG) 
in magnetically trapped atomic vapours 111 has revived 
interest in the theoretical study of Bose gases. Mean-field 
methods provide us with relatively simple predictions 
both for the equilibrium properties of these systems (en- 
ergy per particle, density profiles, condensate fraction) 
and for the dynamic behavior (frequency of collective ex- 
citations, interference effects), which have been found in 
close agreement with experiments (for a review see Ref. 
lyi). In fact, the atomic clouds realized in experiments are 
very dilute, the average distance between particles being 
significantly larger than the range of interatomic forces, 
and mean-field approaches are well suited. However, the 
investigation of effects beyond mean-field theory is an im- 
portant task, which would make these systems even more 
interesting from the point of view of many-body physics. 
Theoretical studies of these effects have already been pro- 
posed, either by analytic inclusion of fluctuations around 
mean-field y,Q or through numerical calculations based 
on quantum Monte Garlo methods H and, more recently, 
also on correlated basis function approaches Q . All these 
investigations are based on the idea that, for the values 
of density relevant in experiments, the details of the in- 
teratomic potential can be neglected and one can safely 
use the hard-sphere model in numerical simulations, and 
the expansion in powers of the gas parameter no?' , fixed 
by the number density n and the s-wave scattering length 
a, in analytic corrections beyond mean-field. The main 
motivation of the present work is to verify the validity of 
this approach. By using a diffusion Monte Carlo (DMC) 
method we calculate the ground-state energy of a system 
of bosons interacting through different two-body model 
potentials. We explicitly show that for the values of the 
gas parameter reached in magnetic traps [no?' ~ lE-5 - 



lE-4) the behavior is universal and fixed by na^ and that 
the corrections to the mean-field energy are much larger 
than the differences due to the details of the interatomic 
potential. 

The ground state of a homogeneous dilute Bose gas was 
intensively studied in the 50 's and early 60 's. One of the 
main results of this investigation is that the ground-state 



energy can be expanded in powers of 



In units of 



h /2ma the energy per particle takes the form 
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The first term corresponds to the mean-field prediction 
and was already calculated by Bogoliubov 0. The cor- 
rections to mean-field have been obtained using pertur- 
bation theory: the coefficient of the (na^)'^/^ term was 
first calculated by Lee, Huang and Yang m, while the 
coefficient of the last term was first obtained by Wu M . 
Both of them were originally derived for hard spheres, 
but it was shown that the same expansion is valid for 
any repulsive potential with scattering length a pO| , pl| . 
Furthermore, Hugenholtz and Pines [ll| have shown that 
the higher-order terms in the expansion (|l|) depend on 
the "shape" of the interatomic potential. It is worth 
pointing out that, for the values of na'^ relevant in exper- 
iments, the corrections to the mean-field energy are very 
small (~ 3%). In a recent paper Lieb and Yngvason |lj] 
have also provided a rigorous lower bound for the ground- 
state energy of a Bose gas holding for non-negative, fi- 
nite range, spherical, two-body potentials. The lower 
bound coincides with the Bogoliubov mean-field term: 
E/N > ATT{na^)h,/2ma^. However, what is not well es- 
tablished up to now is the range of validity of the univer- 
sal law (|l|). 

In realistic systems the above restrictions for the in- 
teratomic potential do not hold, because an attractive 
tail is generally present. In this case the situation for the 
ground state is different. If the potential does not sustain 
any many-body bound state and the scattering length is 
positive, the ground state of the system still behaves like 
a gas and the expansion Im) should hold. Conversely, in 
the case of potentials which have two, three or more-body 
bound states, the ground-state of the system is no longer 
a homogeneous gas, but a state with clusters of atoms 



formed. However, if the scattering length is positive, a 
uniform gas can still exist as a metastable state which 
will be long-lived at very low densities. Our results for a 
simple model of interatomic potential with an attractive 
part show that the energy of the gas- like metastable state 
still follows the universal law given by (0). 

We consider a system of N spinless bosons with 
mass m described by the many-body Hamiltonian H = 
~in^/2m) E.=i V|+E.<, ^(|r-r, I), where V{r) is the 
two-body, spherical, interatomic potential. The DMC al- 
gorithm solves exactly, apart from statistical uncertainty, 
the iV-body Schrodinger equation for the ground-state 
energy of the system (for further details on the method 
see for example Ref. [^). We have used different choices 
for the potential V{r): 
1) Hard-sphere (HS) potential defined by 
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where the diameter a of the hard sphere corresponds to 

the scattering length. 

2) Soft-sphere (SS) potential defined by 



V{r) = 



Vo ir<R) 
{r> R) 



(3) 



with a scattering length: a = R[l — ta,nh{KoR) / KqR] 
with Kq = Vom/h^ and Vo > 0. For finite Vq one has 
always R > a, while for Vq — > +00 the SS potential 
coincides with the HS one with R ~ a. For the SS 
potential we have considered two choices: R = 5a and 
R = 10a. It is worth noticing that the hard sphere and 
the very soft sphere with R = 10a represent two extreme 
cases for a repulsive potential. In the HS case, the en- 
ergy is entirely kinetic, while for the very "soft" potential 
a ~ (m/h ) L V{r)r^dr, according to Born approxima- 
tion, and the energy is almost all potential. In the latter 
case the ground-state wave function is close to the non- 
interacting one (see Fig. 4). 
3) Hard-core square-well (HCSW) potential defined by 



-fcx) (r < Re) 
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for which the s-wave scattering length is given by: a = 
Rc + iR- Rc){l - tan[Ka{R ~ Rc)]/Ko{R - Re)} with 
Ki = Vam/h? and Vq > 0. For fixed Re and R the 
scattering length coincides with the hard-core radius Re 
for Vo = and by increasing Vq it exhibits resonances 
each time a two-body bound state appears in the well. 
In many real potentials (such as in ^^Rb and ^^Na) the 
scattering length is much larger than the size of the atom. 
To reproduce this situation we have chosen R = 5i?c and 
a = lOR with only one two-body bound state in the well. 



In the case of purely repulsive potentials we apply di- 
rectly the DMC algorithm to obtain the ground-state en- 
ergy. As importance sampling we use a Jastrow trial 
function: -0t(R) = -f/Tlri, ■•, tat) = Y{i<j firij)- The 
Jastrow factor /(r) is chosen as the exact wave func- 
tion of a pair of particles interacting through V{r) with 
total energy e for r < R, and for r > R a function 
f{r) = 1 — Ae~^^°' which goes rapidly to one. The pa- 
rameter a is left as a variational parameter, while the 
coefficient A and the matching point R are chosen so 
that /(r) and its first derivative be continous at r = ^ 
and the local energy (— ?i^V^/2m -t- V)f{r)/f{r) be also 
continous at r = _R. The energy e is the second varia- 
tional parameter. Before starting the DMC calculation, 
we perform a variational Monte Carlo (VMC) analysis 
to optimize the parameters of the trial wave function. 
At low densities, we find that the DMC calculation im- 
proves very little on the VMC result. For example, for 
the HS potential at na^ =lE-5 and in units of h /2ma^: 
EvMc/N =1.278(l)E-4 and Edmc/N =1.274(1)E-4. 
This means that at low densities ^px (R-) has a large over- 
lap with the "true" ground-state wave function. 

The case of the HCSW potential needs a careful treat- 
ment. Since the potential has a two-body bound state, 
the gas-like state is not the ground state. To obtain the 
energy of the metastable gas-like state it is necessary 
to project out the bound-state component of the wave 
function. This can be achieved by using the same trial 
wave function as for the SS potential (any trial func- 
tion which is positive in the region where the potential is 
attractive would be equally appropriate) and then pro- 
jecting the results for the energy by means of an auxil- 
iary function '0p(R) ortogonal to any many-body state 
with bound pairs. In the Monte Carlo formulation this 
is realized by the weighted integral E = /(iR/(R, t — > 

oo)[7?Vp(R)]/V'T(R)//dR/(R,r ^ (X3)^p(R)/^t(R), 
where /(R, r — > 00) is the density of walkers generated by 
the DMC calculation. The projecting function V'p(R) is 
chosen as: t/jpCR.) — Yl^^- fpi^ij), where fp{r) coincides 
with the trial two-body function /(r) for r > R, while 
for r < ^ is given by the exact solution for two particles 
interacting through V{r) with energy e > (we require 
the same continuity conditions at the matching point R 
as for f{r)). Since the matching point R is much larger 
than the range R of the potential, fp{r) is orthogonal to 
the bound-state wave function, and V'p(R) is orthogonal 
to any many-body state with one or more bound pairs 
formed. In this way, we eliminate from the calculation all 
states with bound pairs. At low densities we expect that 
bound pairs give the main contribution to the ground 
state compared to three or more-body bound states. By 
eliminating these bound pairs we are thus left with the 
gas-like metastable state with lowest energy. 

We are now in a position to discuss our results. In Ta- 
ble I we present the DMC results for the equation of state 



of the HS potential. In all the calculations we have used 
a simulation box containing TV = 500 particles in order 
to reduce finite size effects well below statistical uncer- 
tainty. Previous Monte Carlo calculations of the ground- 
state energy of the hard-sphere boson system have been 
performed at densities typical of hquid ^He {na^ ~ 0.1) 
M,15[. The results for the two highest densities in Ta- 
ble I are compatible with the Green's Function Monte 
Carlo results of Ref. ||l5| . In Fig. 1 , our results for the 
HS equation of state are shown together with the various 
terms of the analytic expansion (|l|). One can see that the 
Lee-Huang- Yang (LHY) correction [second term in P)] 
represents a significant improvement on the mean-field 
prediction and the inclusion of this term allows for a good 
approximation of the equation of state up to very high 
densities. On the contrary, the logarithmic correction 
[third term in (p] goes wrong already at intermediate 
densities {na^ ~ lE-3). 

In Table II, we report the results of the comparison be- 
tween the various potentials considered in this work (the 
corresponding results for the HS potential can be read 
from Table I). At low values of the gas parameter na^ all 
potentials give the same results and only at the largest 
density reported {na^ =lE-3) the results for the SS po- 
tential start to deviate from the HS values. The univer- 
sal behavior is better shown in Fig. 2, where we plot the 
difference between the calculated energy per particle and 
the mean-field term and compare it with the LHY correc- 
tion. It is worth mentioning that this difference is always 
positive in agreement with the lower bound provided by 
the mean- field term, as discussed in Ref. [Q. The very 
low density regime and the relevance of the logarithmic 
correction is analysed in the inset of Fig. 2. The HS 
and SS results show evidence of the presence of this log- 
arithmic correction, however the effect is tiny and goes 
rapidly wrong for larger densities. Due to larger error 
bars no conclusion can be drawn for the HCSW results 
concerning the logarithmic term. 

Another quantity of interest is the condensate fraction 
No/N which we calculate from the long-range behavior 
of the one-body density matrix: Nq/N — linir^oo pif) 
(see Ref. |l3| for further details). In Fig. 3, we show the 
results for Nq/N as a function of na^ for the HS and the 
two SS potentials using the extrapolated estimator. The 
results are compared with the analytic expansion 
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calculated by Bogoliubov [Q . At low densities this law is 
universal and agrees very well with the results of the three 
model potentials. First deviations from universality start 
to appear for na^ ~ lE-3, result which coincides with the 
emergence of deviations in the energy values. 

In Fig. 4 we show the two-body distribution function 
(;(?■) for hard spheres at various densities obtained from 



the DMC calculation using the method for pure estima- 
tors of Ref. [|l^. One can clearly see the building up 
of short range correlations at high densities. The dis- 
tribution function g{r) can be estimated at low densi- 
ties using the Bogoliubov result for the static form fac- 
tor: S{k) = h'^k^/2me{k), where e(fc) = {h'^k'^/ATn^ + 
gnh^k^/mY^^ is the usual Bogoliubov dispersion rela- 
tion with coupling constant g — Aj^Ti a/m. The Fourier 
transform of S{k) gives directly the distribution function 
g{r), which for r ^ a exhibits the asymptotic behav- 
ior g{r) ~ 1 - [47r^/2(na^)3/2(r/a)^]"^ In Fig. 4 we 
also compare, at the density na^ =lE-4, the distribu- 
tion function obtained from the Bogoliubov S{k) with 
the g{r) of hard spheres and soft spheres with R = 10a. 
The long range behavior is universal and agrees with the 
Bogoliubov result, while at short distances the HS and 
SS potentials give completely different results since in the 
second case particles can penetrate the repulsive poten- 
tial. At higher densities the Bogoliubov approximation 
for g{r) becomes poorer and misses completely the shell 
structure of the distribution function. 

In conclusion, we have shown that for the values of den- 
sity realized in magnetic traps the corrections to mean- 
field are fixed only by the gas parameter na^ and do not 
depend on the details of the interatomic potential. These 
effects, although small, have a simple theoretical descrip- 
tion in terms of the parameter na^ and, possibly, can be 
singled out in future precision measurements H]. 
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TABLE n. Energy per particle for the SS and HCSW po- 
tential (in units of h^ /2ma^) 
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FIG. 1. Equation of state for the HS potential. The solid 
circles are the DMC energies (error bars are smaller than 
the size of the symbols); the solid line corresponds to the 
mean- field prediction [first term in eq. Imj]; the long-dashed 
line includes the LHY correction [first two terms in eq. (dI)]; 
the short-dashed line includes also the logarithmic correction 
and corresponds to the full expansion (fO). 



FIG. 2. Corrections to the mean-field energy. Circles: HS 
potential; down triangles: SS potential {R — 5a); squares: 
SS potential {R = 10a); up triangles: HCSW potential; solid 
line: LHY correction [second term in eq. (Fy)]. Error bars 
are smaller than the size of the symbols. The inset shows the 
results for the HS, SS {R = 10a), and HCSW potentials in 
the extremely low density region, the solid line is again the 
LHY correction, while the dashed line includes the logarith- 
mic correction [second and third term in eq. m)]. 



FIG. 3. Condensate fraction as a function of the gas pa- 
rameter. Circles: HS potential; down triangles: SS potential 
{R = 5a); squares: SS potential {R = lOo). The solid line 
corresponds to the Bogoliubov expansion (H). Error bars are 
smaller than the size of the symbols. 



FIG. 4. Two-body distribution function. Solid lines: 
hard-spheres at densities na'^ =lE-4 (lowermost), lE-2, 0.1, 
0.244 (uppermost). Long-dashed line: soft-spheres {R — 10a) 
and short-dashed line: Fourier transform of the Bogoliubov 
static form factor S{k) at the lowest density na^ =lE-4. 



TABLE I. Energy per particle for the HS potential 
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